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A novel sampling theorem on the sphere 
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Abstract — We develop a novel sampling theorem on the sphere 
and corresponding fast algorithms by associating the sphere 
with the torus through a periodic extension. The fundamental 
property of any sampling theorem is the number of samples 
required to represent a band-limited signal. To represent exactly 
a signal on the sphere band-limited at L, all sampling theorems 
on the sphere require 0{L'^) samples. However, our sampling 
theorem requires less than half the number of samples of 
other equiangular sampling theorems on the sphere and an 
asymptotically identical, but smaller, number of samples than 
the Gauss-Legendre sampling theorem. The complexity of our 
algorithms scale as 0{L ), however, the continual use of fast 
Fourier transforms reduces the constant prefactor associated with 
the asymptotic scaling considerably, resulting in algorithms that 
are fast. Furthermore, we do not require any precomputation 
and our algorithms apply to both scalar and spin functions 
on the sphere without any change in computational complexity 
or computation time. We make our implementation of these 
algorithms available publicly and perform numerical experiments 
demonstrating their speed and accuracy up to very high band- 
limits. Finally, we highlight the advantages of our sampling 
theorem in the context of potential applications, notably in the 
field of compressive sampling. 

Index Terms — Harmonic analysis, sampling methods, spheres. 



I. Introduction 

IN many fields of science and engineering data are measured 
on a spherical manifold. Applications where data are de- 
fined inherently on the sphere are found in computer graphics 
(e.g. [1]), planetary science (e.g. [2]-[5]), geophysics (e.g. 
[6]-[8]), quantum chemistry (e.g. [9], [10]) and astrophysics 
(e.g. [11], [12]), to quote only a few. In many of these 
applications a harmonic analysis of the data is insightful. For 
example, spherical harmonic analyses have been remarkably 
successful in cosmology over the past decade, leading to the 
emergence of a standard cosmological model. Observations of 
the anisotropics of the cosmic microwave background (CMB), 
which are made on the celestial sphere, contain a wealth of 
information about the early Universe. Cosmologists extract this 
information from the angular power spectrum of observations 
of the CMB, computed through a harmonic transform on the 
sphere (e.g. [13]). Recent and upcoming full-sky observations 
of the CMB are of considerable size, containing approximately 
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three [12] and fifty [14] million samples respectively. Fur- 
thermore, observations are made of both the temperature and 
polarisation of the CMB, which give rise to scalar and spin ±2 
functions on the sphere respectively. Consequently, the ability 
to perform fast scalar and spin spherical harmonic transforms 
on large data sets is of considerable importance in cosmology 
and beyond. 

Algorithms to perform spherical harmonic transforms have 
received considerable attention already. Some correspond to 
sampling theorems on the sphere, where the forward and 
inverse transform are theoretically exact for a band-limited 
signal on the sphere. Others adopt approximate quadrature 
rules on the sphere, resulting in approximate harmonic trans- 
forms that do not correspond to sampling theorems on the 
sphere. However, these approximate algorithms typically arise 
from particular pixelisations of the sphere which meet de- 
sirable practical criteria, such as pixels of equal area, and 
are no less important. We focus here on approaches that 
lead to sampling theorems on the sphere with theoretically 
exact transforms for signals on the sphere band-limited at 
L (where the harmonic band-limit L is defined formally in 
Sec. |III-B| l. Note that the cuiTent [12] and forthcoming [14] 
CMB observations discussed previously support band-limits 
of L = 1024 and L = 4096 respectively. All sampling 
theorems require N samples on the sphere of order 0{L'^), 
however the exact number of samples required varies for each 
sampling theorem. For many applications reducing the number 
of samples required to represent a band-limited signal on the 
sphere is of fundamental importance. 

Sampling theorems on the sphere and their associated 
numerical algorithms are evaluated by four criteria: (i) the 
number of samples required to represent a band-limited signal 
exactly; (ii) their computational complexity; (iii) their speed; 
and (iv) issues surrounding any precomputation. From an 
information theoretic viewpoint, the fundamental property of 
any sampling theorem is the number of samples required 
to represent a band-limited signal exactly. In this article we 
review algorithms to compute spherical harmonic transforms 
accurately and efficiently. We also present a novel sampling 
theorem on the sphere and corresponding fast algorithms. 
Our approach compares favourably to the state-of-the-art as 
evaluated by the four criteria listed previously. Furthermore, 
our algorithms apply to both scalar and spin functions on 
the sphere without any change in asymptotic complexity or 
computation time. 

The remainder of this article is structured as follows. In 
Sec. [n] we review comprehensively the literature regarding 
the computation of spherical harmonic transforms, placing our 
novel sampling theorem and fast algorithms in the context of 
preceding work. Harmonic analysis on the sphere is reviewed 



concisely in Sec. Ill to present the mathematical preliminaries 
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required subsequently. Our sampling theorem on the sphere 
and the corresponding fast algorithms to compute spherical 
harmonic transforms are derived in Sec. |IV] In Sec. |V] we 
evaluate our algorithms numerically and discuss the advan- 
tages of our sampling theorem in the context of potential 
applications, notably in the field of compressive sampling. 
Concluding remarks are made in Sec. |VT] 

II. Review 

The development of sampling theorems on the sphere and 
fast algorithms to compute spherical harmonic transforms has 
been driven largely by researchers in the fields of computa- 
tional harmonic analysis, geophysics and astrophysics. Due 
to the diverse nature of these fields, the literature on this 
topic appears to be somewhat disjoint. We attempt to unify 
these works here and to present a comprehensive review of 
the historical development of the field. 

For isolatitude sampling schemes, where the samples are 
gathered in isolatitude annuli, a separation of variables may be 
performed to rewrite the scalar spherical harmonic transform 
as a Fourier transform in longitude and an associated Legendre 
transform in colatitude. The computation is then dominated by 
the associated Legendre transform, reducing the complexity 
from ©(L^) to 0{L^). Isolatitude samplings are thus very 
common and hence we restrict our attention to such schemes. 
Spin lowering and raising operators can be used to relate 
the harmonic transform of spin functions to the scalar case, 
hence scalar transforms have received the majority of attention. 
Approaches to improve the computational performance of 
scalar spherical harmonic transforms attempt either to reduce 
the asymptotic complexity further through fast associated Leg- 
endre transforms or to reduce the constant prefactor associated 
with the asymptotic scaling of the algorithm. 

First attempts to compute a scalar spherical harmonic trans- 
form through a fast Legendre transform were performed by 
Orszag [15] and were based on a Wentzel-Kramers-Brillouin 
(WKB) approximation. Alternative approaches using the fast 
multipole method (FMM) [16] have been considered by Alpert 
& Rokhin [17] and by Suda & Takani [18]. The complexity 
of these algorithms scale linearly with the desired accuracy. 
An alternative approximate algorithm using WKB frequency 
estimates has been developed by Mohlenkamp [19], [20] for 
functions with band-limits that are a power of two, however 
this approximation can be controlled independently of com- 
plexity, which scales as 0{L^ log2-Zj). In any case, these types 
of approach are necessarily approximate and do not yield exact 
sampling theorems on the sphere. 

Other approximate approaches based solely on the sepa- 
ration of variables have been developed more recently for 
pixelisations of the sphere that satisfy certain practical require- 
ments, such as HEALPi>|^ [21] and IGLOC^[22], resulting 
in algorithms of complexity O(L^). For these pixelisations 
only approximate quadrature rules exist, hence the spherical 
harmonic transform algorithms of HEALPix and IGLOO 
are not theoretically exact. Nevertheless, these pixelisation 



schemes satisfy a number of desirable practical criteria, such 
as pixels of equal area, and their associated harmonic trans- 
form algorithms are of sufficient accuracy for many practical 
purposes. These schemes have found considerable application 
in the analysis of CMB data. 

Exact transforms with associated sampling theorems have 
been constructed for particular pixelisation schemes. It is 
well-known that Gauss-Legendre quadrature may be used to 
construct exact spherical harmonic transforms. To our knowl- 
edge, this result was first highlighted in published work by 
Shukowsky [23], which in turn refers to unpublished (and in- 
accessible) work by Payne from 1971 [24]. An exact sampling 
theorem can be constructed from A^gl — L{2L — 1) ~' 2L^ 
samples on the sphere, where the sample locations in colati- 
tude are chosen as the roots of the Legendre polynomials of 
order L, as dictated by Gauss-Legendre quadrature. Through a 
separation of variables, the resulting algorithm is 0{L^). The 
GLESl|^[25] pixelisation scheme has been constructed using 
Gauss-Legendre quadrature, however this scheme uses twice 
as many samples in colatitude as required, i.e. approximately 
2A^GL ^ samples are used. GLESP has also found 

considerable application in the analysis of CMB data. 

The first theoretically exact sampling theorem on an equian- 
gular pixelisation was developed by Shukowsky [23], requiring 
A^s — {"^L — 1)^ ~ 4L^ samples on the sphere; however, the 
exactness of this approach was not studied numerically. Al- 
though this algorithm remains 0{L'^), a separation of variables 
may be used to reduce the computational complexity to 
0{L^). An alternative sampling theorem on the sphere for an 
equiangular pixelisation was developed by Driscoll and Healy 
[26]. Moreover, a divide-and-conquer approach to computing 
a fast associated Legendre transform in the cosine basis was 
derived [26]. The resulting algorithm is exact in exact precision 
arithmetic, and has computational complexity ©(L^logli), 
but is known to suffer from stability problems [27], [28]. 
Healy et al. [27] readdressed the work of Driscoll & Healy 
[26], reformulating the sampling theorem on the sphere and 
developing some variants of the original algorithm]^ which 
are available for download]^ However, the only variant that is 
universally stable is the so-called semi-naive algorithm, which 
remains 0{L^). Algorithms to compute spin ±2 transforms 
are derived by Wiaux et al. [29] and Kostelec et al. [30] using 
spin raising and lowering operators to relate spin transforms 



http : / /healpix . jpl . nasa . gov/ 




^ http : //www .mrao . cam. ac . uk/pr 


ojects/cpac/ igloo/ 



^http: / /www.glesp.nbi .dk/| 
Healy et al. [27] derive a number of variants of the original Driscoll 
& Healy algorithm [26], including the so-called semi-naive, .simple-split and 
hybrid algorithms. The semi-naive algorithm avoids dividing (and conquering) 
the problem, resulting in 0{L'^) complexity. The simple-split algorithm is 
a simpler and more stable divide-and-conquer approach than the original 
algorithm but with an increased complexity of C>{L^/^ \ofr}J^L) and is less 
stable than the semi-naive approach. The splitting required by the simple- 
split algorithm is costly (in terms of execution time rather than asymptotic 
complexity), thus for a band-limit of L = 1024 the semi-naive algorithm 
is greater than two times faster than the simple-split algorithm [27]. The 
hybrid algorithm attempts to mitigate the slow execution of the simple-spit 
algorithm and the higher complexity of the semi-naive algorithm by splitting 
the problem between them. The hybrid algorithm appears to achieve a good 
compromise between stability and efficiency. However, the overall complexity 
of this algorithm is not clear since it depends on the split between the semi- 
naive and simple-split algorithms and on user sp ecified pa rameters. 

^http : //www. cs ■ dartmouth.edu/-geelong/ sphere/] 
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to the scalar case, before applying the scalar algorithms de- 
veloped by Healy et al. [27]. In general, this type of approach 
may be used to compute spin transforms for arbitrary spin 
[31], however the complexity of the resulting algorithm then 
also scales linearly with spin number All of these algorithms 
[26], [27], [29], [31], [32] require A^dh = 2L(2L - 1) ~ 
samples on the sphere and, moreover, the divide-and-conquer 
based approaches are all restricted to harmonic band-limits that 
are a power of two. Furthermore, all of the methods with com- 
plexity below O(L^) require a precomputation that requires 
0{L'^) computations and storage. At band-limit L = 1024, for 
example, the precomputation requires 1.2GB of storage [29], 
scaling to approximately 77GB for the band-limit L = 4096 
of forthcoming CMB observations. Precomputation quickly 
becomes infeasible for high band-limits, thus the O(L^) semi- 
naive algorithm is the most universally applicable fast algo- 
rithm implementing the DriscoU & Healy sampling theorem. 

Other approaches to reduce the cost of computing spherical 
harmonic transforms focus on reducing the constant prefac- 
tor associated with the asymptotic complexity of algorithms. 
These approaches have typically exploited fast Fourier trans- 
forms (FFTs) on equiangular pixelisations to reduce compu- 
tation time through an association between the sphere and the 
torus, while their complexity remains 0{L^). To our knowl- 
edge, the first algorithm based on this technique was developed 
by Dilts [33], where the North and South poles of the sphere 
were identified to map the sphere to the torus. However, this 
algorithm is approximate and does not result in a sampling 
theorem on the sphere. One of the authors of this article 
developed a samphng theorem on the sphere [34] by making 
periodic extensions of the sphere in colatitude in order to make 
an association with the torus. This sampling theorem requires 
A^M = - 1)(2L - 1) + 1 samples on the sphere 
and, moreover, applies to any spin number without the applica- 
tion of spin lowering and raising operators. However, the for- 
ward algorithm associated with this samphng theorem proved 
imstable (the inverse algorithm did not suffer from stability 
issues). Recently, Huffenberger & Wandelt [35] adopted the 
inverse transform of this approach and resolved the instability 
in the forward algorithm, although in doing so increased the 
number of points required to sample a band-limited function 
on the sphere to A/hw = 2L(2L — 1) ~ 4i^. In this article we 
readdress sampling theorems derived by associating the sphere 
with the torus through periodic extensions and develop a sam- 
pling theorem requiring A\iw = (L— 1)(2L — 1) + 1~ 2L? 
samples on an equiangular pixelisation, with corresponding 
fast algorithms that do not suffer from any stabihty issues. 



III. Harmonic Analysis on the Sphere 

In this section we review harmonic analysis on the two- 
sphere S^. We first review the scalar spherical harmonic 
transform, before generalising to the spin case. Associations 
are then made between the spin spherical harmonics and the 
Wigner functions, where the latter provide an orthogonal basis 
for the decomposition of square integrable functions on the 
rotation group SO (3). 



A. Scalar spherical harmonics 

We consider the space of square integrable functions on the 
sphere L^(S^), with the inner product of f,g € L^(S^) defined 

{f,g)= / dm^)mip)g*{e,ip), 

where dfl{9, ip) = smOAOAip is the usual invariant measure 
on the sphere and (6*, ip) define spherical coordinates with 
colatitude 9 € [0,7r] and longitude € [0,27r). Complex 
conjugation is denoted by the superscript *. 

The scalar spherical harmonic functions form the canonical 
orthogonal basis for the space of L^(S^) scalar functions on 
the sphere and are defined by 

for natural f e N and integer m & "L, \m\ < I, where 
Pp{x) are the associated Legendre functions. We adopt the 
Condon-Shortley phase convention, with the (—1)™ phase 
factor included in the definition of the associated Legendre 
functions, to ensure that the conjugate symmetry relation 
^/m(^.¥') = {-ir Yc-miO,^) holds. The orthogonality 
and completeness relations for the spherical harmonics read 
{Yem, Yfrn') = Su'Smm' and 

cc £ 

E ^irnio, ^) Y;^{e', = 6{cose- cose') - 

respectively, where Sij is the Kronecker delta symbol and S{x) 
is the Dirac delta function. 

Due to the orthogonality and completeness of the scalar 
spherical harmonics, any square integrable scalar function on 
the sphere / e L^(S^) may be represented by its spherical 
harmonic expansion 

oo i 
£=Q m=-i 

where the spherical harmonic coefficients are given by the 
usual projection onto each basis function: f^m — (/, X^m)- 
The conjugate symmetry relation of the spherical harmonic 
coefficients of a real function is given by = (— 1)™ fe,-m, 
which follows directly from the conjugate symmetry of the 
scalar spherical harmonic functions. 

B. Spin spherical harmonics 

Square integrable spin functions on the sphere sf S L^(S^), 
with integer spin ,s G Z, are defined by their behaviour under 
local rotations. By definition, a spin function transforms as 

,/'(0,^)=c--^,/((?,^) (1) 

under a local rotation by x, where the prime denotes the 
rotated function. It is important to note that the rotation 
considered here is not a global rotation on the sphere, such as 
that represented by an element of the rotation group S0(3), 
but rather a rotation by x in the tangent plane at {6, ip). The 
sign convention that we adopt here for the argument of the 
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complex exponential in ([T]) differs to the original definition 
[36] but is identical to the convention used recently in the 
context of the polarisation of the CMB [37]. 

The spin spherical harmonics sYim{&,f) form an orthog- 
onal basis for L^(S^) spin s functions on the sphere for 
|s| < £. Spin spherical harmonics were first developed by 
Newman & Penrose [36] and were soon realised by Gold- 
berg [38] to be closely related to the Wigner functions. We 
therefore defer the explicit definition of the spin spherical 
harmonic functions until Sec. |III-C| The conjugate symme- 
try relation given for the spin spherical harmonics is given 
by sY;^id,ip) = (-l)*+™_,r,,_™(0,</j). The spin spherical 
harmonics satisfy identical orthogonality and completeness 
relations as the scalar spherical harmonics. 

Due to the orthogonality and completeness of the spin 
spherical harmonics, any square integrable spin function on 
the sphere sf £ L^(S^) may be represented by its spherical 
harmonic expansion 

oo e 

i=0 m=-l 

where the spin spherical harmonic coefficients are given by the 
usual projection onto each basis function: sfim — (sf, s^^m)- 
Note that the spin spherical harmonics and transforms simply 
generalise the scalar equivalents to spin signals, reducing to 
the standard scalar case for s = 0. When deriving our novel 
sampling theorem we consider signals on the sphere band- 
limited at L, that is signals such that sftm = 0, V£ > L. The 
conjugate symmetry relation of the spin spherical harmonic 



s f^^ — m 



for a 



coefficients is given by s.f*em = (^1 
function satisfying sf* = -sf (which for a spin s ~ 
function equates to the usual reality condition) and follows 
directly from the conjugate symmetry of the spin spherical 
harmonics. 

Spin raising and lowering operators, 3 and § respectively, 
exist so that spin s ± 1 functions may be obtained from spin 
s functions [36], [38]. Spin raising and lowering operators 
are often used repeatedly to relate spin s functions to scalar 
functions on the sphere. 



C. Wigner functions 

The Wigner functions Z?f„„(Q!, /3, 7), for natural £ E N 
and integer m,n E Z, form an orthogonal basis for the 
space L^(S0(3)) of square integrable functions on the rotation 
group, and are parameterised by the Euler angles (a,/?, 7), 
where a E [0, 27r), j3 E [0,7r] and 7 E [0,27r)QThe Wigner 
functions may be decomposed as [39] 



-ima ^£ 



m7 



(2) 



^We adopt the zyz Euler convention corresponding to tlie rotation of a 
pliysical body in a fixed co-ordinate system about the z, y and z axes by 7. 
/3 and a respectively. 



where the real polar d-functions are defined by [39] 

/ {£ + ny.{£~n)l ( . /3 

/ o\ n+m 

x(^cos|j p(:r"^"^"^(cos/3), (3) 

where Pj,'^'''\-) are the Jacobi polynomials. Note that recur- 
sion formulae are available to compute rapidly the Wigner 
d-functions (e.g. [40], [41]). The d-functions satisfy a number 
of symmetry relations; in this work we make use of the 
symmetry relations [39] 



and 



rfinM) = (-l)'"""rfL.(/3)- 



(4) 
(5) 

(6) 



We are not concerned with decompositions of functions on 
the rotation group in this article but rather representations 
of the spherical harmonics by Wigner functions. The spin 
spherical harmonics may be defined by the Wigner functions 
through [38] 



2£+l 
An 



(7) 



Defining the spherical harmonics in this manner allows us to 
apply standard Wigner function decompositions to the spher- 
ical harmonic functions. We subsequently make considerable 
use of the Fourier series decomposition of the d-functions 
given by [42]: 



E 



A , A , e 



(8) 



where A^„ = (i^„(7r/2). This expression follows from a 
factoring of rotations as highlighted by Risbo [40]. The Fourier 
series representation of dm„(/3) given by ([sj allows one to 
write the spherical harmonic expansion of sf in terms of 
a Fourier series expansion of sf extended appropriately to 
the two-torus (as discussed in more detail in Sec. IV 1. 



Consequently, ([8]) is fundamental to the derivation of our 
sampling theorem on the sphere and fast algorithms. 

IV. Fast Spherical Harmonic Transform 

We derive fast algorithms for performing forward and 
inverse spin spherical harmonic transforms and discuss the 
corresponding sampling theorem on the sphere. Our approach 
involves an extension of the sphere to the torus so that 
FFTs may be exploited to reduce the cost of computation. 
It is related closely to the algorithms derived by one of 
the authors in a previous work [34] and to the algorithms 
derived by Huffenberger & Wandelt [35], however it does 
not suffer from the instabilities of the former approach and 
requires only half as many samples on the sphere as the latter 
approach. We first present the general harmonic formulation 
of our algorithms, followed by a discussion of periodisation 
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and discretisation, before deriving algorithms to perform exact 
forward and inverse transforms. Our sampling theorem also 
leads to a new quadrature rule on the sphere, which we then 
present, followed by a discussion of symmetries that may be 
exploited to improve the efficiency of computation for real 
signals. 



A. Harmonic formulation 

We consider the harmonic transform of spin functions on 
the sphere sf G L^(S^), band-limited at L\ consequently, all 
summations over or up to £ are truncated to L— 1. Furthermore, 
harmonic coefficients are not defined for \m\ > £, hence we 
define them to be zero to enforce the contraint |m| < £ when 
summations are interchanged. 

By noting the definition of the spin spherical harmonics in 
terms of Wigner functions (|7]|, the Wigner decomposition (|2]l 
and the Fourier expansion of the Wigner d-functions ([8]l, the 
forward transform of s/ may be written 



Jim = (-l)^i' 



21 + 1 



where 



and 



mm' 



m'=-(L-l) 



(9) 
(10) 

(11) 



In Sec. IV-D we consider implicit quadrature rules to evaluate 
( [T0| and exactly. By noting the same substitutions and 
interchanging the order of summation, the inverse transform 
of sfem tie written 



L-l 
m=-(L-l) 



where 



and 



L-l 



s-^m(^) ^ ^ s-^mm' 

m' = -{L-l) 



(12) 



(13) 



L-l 



= (—1)* 



f=0 



2£+l 
An 



m'ni r, 



(14) 

Although recasting the forward and inverse spherical harmonic 
transforms in this manner is no more efficient than the original 



formulation, expressions (10i-(13i highlight similarities with 



Fourier series representations. However, the Fourier series 
expansion is only defined for periodic functions; thus, to recast 
these expressions in a form amenable to the application of 
Fourier transforms we must make a periodic extension in 
colatitude 9. 



B. Periodic extension 

We make a periodic extension of to the domain [0, 27r) 
so that FFTs may be used to compute the forward and 
inverse spherical harmonic transform rapidly. When making 
this periodic extension we must be careful to ensure that 
the symmetry of our current representation is respected on 
the new domain; we must apply the symmetries dictated 
by the inverse transform when imposing periodisation in the 



forward transform. By substituting (12i into (111 and noting 



the continuous orthogonality of the complex exponentials, we 
find the forward and inverse expressions in 9 are related by 



sG,n{e)=2TT sF„,{9) 



(15) 



Consequently, the symmetry we impose in sGm{S) when 
extended periodically must match the symmetry of sFm{d)- 
By reflecting 9, we obtain the following symmetry for sFm{0)'- 

e 

sFfn{ 9^ ^ ^ sFmm' ^ ' ( 1) sFmiP^ 5 



where we have noted the symmetry 



F 

m. — m 



(-1)™+" F , 

\ mm 



(16) 



following from (14 1 and (jsj. Thus, we extend sGm{d) to the 
[0, 27r) domain by constructing^ 



sGm{Q) — 



sGm{9) , 

(-l)™+%G,n(2^ 



9 e [0, tt] 
9 e (7r,27r) 



Note that we adopt a different periodic extension to other 
approaches framed on the torus [34], [35] by applying the 
extension to the Fourier transform of sf in tp, i.e. to sGm{&)- 
Two periodic extensions, one even and one odd, were required 
in the approach taken previously by one of the authors [34]. 
Huffenberger & Wandelt [35] apply the (—1)™ factor as a shift 
in If by tt, removing the need for two periodic extensions 
but requiring a even number of samples in cp, precluding 
an association with the odd number of points in m unless 
oversampling is performed. We avoid these restrictions by 
applying the periodic extension to the Fourier transform in 
(f of sf, rather than to sf directly. 

C. Discretisation 

We adopt an equiangular sampling of the sphere with 
sample positions given by 

a _ 7r{2t+l) 



2L- 1 



where t e {0, 1, . . . , L - 1} 



(17) 



'We check that this periodic extension does not impose disconti- 
nuities at the poles 9* £ {0,7r}. To avoid discontinuities we require 
sFm{djL) = {—i)"^^" sFm{d*), which follows trivially for m + s even. 

) = _^(6»*), which, due to the 

for m -|- s 

odd. Combined with the representation, for m + s odd. 



From |6| we find _^ 

continuity of the Wigner d-functions, implies df^ _g{9*) 



/ 2/ -I- 1 

this implies sFm{6*) = for m + s odd; hence our periodic extension does 
not impose any discontinuity. 
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and 



2ttp 
2L-1'- 



where {0, 1, . . . , 2L - 2} 



(18) 



In order to extend the 6 domain to [0, 2tt) we simply extend 
the domain of the t index to include {L, L + 1, . . . , 2L — 1}. 

An odd number of sample points are required in both 9 and 
(fi in the extended domain so that a direct association may 
be made with the harmonic indices m and to', resulting in 
A^MW = [L — l)(2i — 1) + 1 ^ 2L? samples on the sphere. 
We also require a symmetric sampling in about the South 
pole so that samples on the extended domain can be obtained 
by reflecting samples defined on the original domain. The node 



positions specified by (17 1 and (18 i eliminate repeated samples 
at the poles 9 — and 9 — 2ti since these points are excluded 
from the pixelisation. However, it is not possible to eliminate 
repeated samples at 6 = tt, since we require a discretisation 
that is symmetric about tt but which contains an odd number 
of sample points. 



D. Forward transform 

The algorithm we derive in this section to compute a 
forward spin spherical harmonic transform essentially follows 
the harmonic formulation presented in Sec. |IV-A| however we 
discuss implicit quadrature rules for the exact evaluation of 
([T0| and ([TT). 

Since ([TTJ is simply a Fourier transform we may appeal 
to the discrete and continuous orthogonality of the complex 
exponential to express this integral exactly by 



27r 



2L- 1 



L-l 

E 

--(L-l) 



for t e {0, 1, . . . , i — 1}. An FFT may be used to compute 
sGm{Ot) for all TO and t with computational complexity 
©(L^logsi). We extend sGm{9) to the domain 9 e [0, 27r) 
through the construction 



i G {0,1,...,L-1} 
■2L-2-t), te{L,...,2L-2} ' 

noting 271- 9t = 62L-2-t- 

We now consider an implicit quadrature rule for the exact 




evaluation of sGmm' through ( 10 1. Firstly, however, we com- 
pute sFmm' from sG„i{9t) by noting \l5\ and by appealing 
to the discrete orthogonality of the complex exponentials to 
invert ([TSj, giving 



L-l 



27r(2L-l) ^ 

^ ' t=-{L-l) 



sGr, 



-'un'Ot 



An FFT may be used to compute s-F^m' for ^ ^nd m' 
with computational complexity 0{L'^ log2 L). Now that we 



to yield 



mm' 



AO iiin9 sGyn{0) e~ 



27r ^ sFmm" w{m" - to') , 

m" = -(L-l) 



(19) 



where the weights are given by 

w{m')= f d9sm9e"'' 
Jq 




to' = ±1 

to' odd, m' 7^ ±1 



2/ (1 — to' ), to' even 



Note that the definition of these weights is identical to that 
derived by Huffenberger & Wandelt [35], however we correct 
some (typographical) errors in their explicit evaluation. Since 
we are concerned with the values of sGmm' for |m'| < L — 1, 
the computation of sGmm' through ( [T9] l explicitly requires 
weights with argument up to ±2(L — 1). If the range of to' 
is extended to |m'| < 3(i — 1), (19 1 may be seen as a (re- 



have sFmm' to hand, we substitute ( 13 i into ( [T0| , noting ( 15 1, 



fleeted) convolution, which may be computed more efficiently 
following a Fourier transform. However, since only the range 
|m'| < L — 1 is of interest, aliasing may be tolerated provided 
that it is outside of this range. To ensure that this is the case 
we zero-pad sFmm' in the domain |to'| e {L, . . . ,2L — 2} 
prior to computing an inverse Fourier transform. We then 
compute an inverse Fourier transform of the weights on the 
same extended domain and take the product of these terms, 
the Fourier transform of which gives sGmm'- Using FFTs to 
compute sGmm' in this manner reduces the computational 
complexity from 0{L^) for the direct calculation of (19i to 
0(L2log2 L). 

Once we have computed sGmm' through the implicit 
quadrature rule discussed previously, we simply compute the 
spin spherical harmonic coefficients sfem through (|9|. The 
complexity of this computation is 0{L^), which dominates 
the overall complexity of the forward algorithm. 

The forward algorithm may be summarised conceptually as 
follows, where we view the upsampling and application of 
weights in the spatial domain: 
1: procedure Forward Transform( gf ) 
2: compute the Fourier transform of s/ in 
3: extend the resultant function to 27r in 9 
4: upsample the resultant function in 6 
5: multiply by the inverse Fourier transform of the 
reflected weights and take the Fourier transform 
in 9 to give the coefficients sGmm' 
6: compute the spherical harmonic coefficients sf im 

from sGmm' 

1: return 
8: end procedure 

Although the complexity of this approach remains identical 
to a standard separation of variables, the continual use of FFTs 
reduces the constant prefactor associated with the asymptotic 
scaling considerably, resulting in an algorithm that may be 
used to compute harmonic coefficients rapidly. Furthermore, 
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a precomputation is not required, which can otherwise ne- 
cessitate very large storage requirements for high band-limits. 
Finally, note that this algorithm applies to both scalar and spin 
functions without any change in computational complexity or 
computation time. 

E. Inverse transform 

The inverse algorithm follows the harmonic formulation 
presented in Sec. IV-A| closely. Firstly, sFmm' is computed 
directly through ([14]), where we also exploit the symmetry 
(16 1 to reduce the number of computations by a factor of 



two. The complexity of this computation is 0{Lr'), which 
dominates the overall complexity of the inverse algorithm. 
FFTs are then used to evaluate ( [T3| and ( [T2| rapidly over 
the extended domain, with complexities 0{Lr\og2 L) for both 
computations. To facilitate the efficient calculation of ( [T2] i 
through the use of an FFT, we compute function values on the 
extended 6 domain [0, 2tt), however we discard those values 
computed in the domain (tt, 27r). The algorithm presented 
here to compute the inverse transform is identical to that first 
proposed by one of the authors [34] and subsequently adopted 
by Huffenberger & Wandelt [35]. The inverse algorithm may 
be summarised as follows: 
1: procedure Inverse Transform( ^/^^ ) 
2: compute the Fourier coefficients sFmm' from sfim 
3: compute the function samples on the extended 

domain by an inverse Fourier transform 
4: construct s/ by discarding samples computed in the 

9 domain (tt, 2tt) 
5: return sf 
6: end procedure 

Since we construct algorithms to perform forward and 
inverse spherical harmonic transforms that are theoretically ex- 
act, our construction corresponds to a novel sampling theorem 
on the sphere. Moreover, to represent a band-limited signal on 
the sphere our sampling theorem requires less than half the 
number of number of samples required by other equiangular 
sampling theorems on the sphere [23], [26], [35] and an 
asymptotically identical, but smaller, number of samples that 
the Gauss-Legendre sampling theorem. 




(a) Weights for L = 4 





(b) Difference for L = 4 





(c) Weights for L = 64 



(d) Difference for L = 64 



Fig. 1. Exact quadrature weights corresponding to our sampling theorem. In 
the left column of panels the weights v{6t) (red squares) defined on [0, 2it) 
and the quadrature weights q{Ot) (yellow diamonds) defined on [0,7r] are 
plotted. These values are compared to samples of the function defined by 
sin(6') on [0, tt) and zero on [tt, 2it) (solid black line). In the right column 
of panels the difference between the quadrature weights q{Ot) and sin(9) are 
plotted. 



may be tolerated, hence the number of samples required in ip 
is reduced from 2L — 1 to L. From the (reflected) convolution 
(19 1, it is apparent that the coefficients sFmm' are required 



for all m' (for m = only). Consequently, the sampling in 9 
remains unchanged, with L samples. However, only weights 
with argument up to ±(X — 1) are required. The (reflected) 
convolution thus spans the range |m'| < 2{L ~ 1). However, 
since only to' = is of interest aliasing may be tolerated 
in m' in all Fourier coefficients sGmm' except to' = 0, 
so that zero-padding is not required before computing the 
(reflected) convolution as a product in the spatial domain. The 



absence of upsampling leads to the explicit quadrature (20i, 
with L(L — 1) + 1 samples on the sphere, where the weights 
are defined by 



F. Quadrature 

The construction of our sampling theorem on the sphere can 
be used to define an explicit quadrature rule for the integration 
of a function band-limited at L. This integration, which 
corresponds to computing the spherical harmonic coefficient 
s/qq, requires approximately half as many samples as needed 
to compute all spherical harmonic coefficients. We define the 
explicit quadrature weights q{dt) to evaluate the following 
integral exactly by the finite sum: 



/ = 



dn{9,^)j{9,^) = 



L-1 L-1 

EE 

t=0 p=0 



where (p'^ = 2ttp/L. As seen from (N^ 
s/qo requires sGmm' for m = to' = I 
aliasing in to in all Fourier coefficients 



J{9uip;)q{9t), (20) 

, the computation of 
) only. Consequently, 
■.Gmm' except m — Q 



2-K 



TlL-2-t) 



and where v{9t) is the inverse discrete Fourier transform of 
the reflected weights w{—m'): 



2L 



1 

E w{-m') 



,Tl' = -(L-l) 



The weights v{^9i) defined on [0, 27r) are exactly the samples of 
the function defined by sin(6') on [0,7r) and zero on [7r,27r), 
band-limited at L. The quadrature weights q{9t) defined on 
[0, tt] are constructed simply by folding the contributions of 
v{9t) on (tt, 2tt) back onto the [0, tt] domain. Both of these 
weights are plotted and compared to sin(0) in Fig. [T] 
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G. Real signals 

In many practical applications signals observed on the 
sphere satisfy a reality condition. For spin signals, the re- 
ality condition is given by gf* = ^sf, which implies the 
conjugate symmetry condition s/|,„ = (—1)"+™ ^sfi.-m on 
the spherical harmonic coefficients of the signal. When this 
reality condition is satisfied (for example, when considering 
the polarisation of the CMB), we may exploit this symmetry to 
recover the harmonic coefficients of a spin ~s signal for free 
from the coefficients of a spin s signal. For the spin s = case, 
the reality condition reduces to the standard reality condition 
of a scalar signal, implying /^^^ = (—1)™ fe,-m- In this case. 



250 p 
200 



noting ( 14 1 and (ml, we obtain the symmetry 



— oF, 



0-^ mm' 



We exploit these symmetries to reduce the computational cost 
of both the forward and inverse algorithms by an additional 
factor of approximately two for real spin s = signals. 

V. Evaluation 

We have implemented our fast algorithms to compute 
spherical harmonic transforms in double precision arithmetic, 
exploiting all of the symmetries discussed in Sec. IV to 
optimise the implementation]^ The core implementation used 
for the numerical experiments presented in this section is 
written in C, using the FFT1/\|^ package to compute Fourier 
transforms, however we also provide a MAT LAB interface. 
We make our Spin Spherical Harmonic Transform (SSHT) 
package containing this implementation available publiclyf*"] 
For comparison purposes, we also implemented in the SSHT 
package an optimised algorithm to compute spherical har- 
monic transforms for the Gauss-Legendre sampling theorem 
on the sphere]^ This algorithm is based on a separation 
of variables and a direct application of the Gauss-Legendre 
quadrature rule, resulting in complexity 0{L^). 

In this section we evaluate our sampling theorem and fast 
algorithms in terms of the number of samples required to 

*The structure of our algorithms suggest multiple transforms of different 
spin may in theory be computed simultaneously at lower cost than consecutive 
computation (since Wigner d-functions do not need to be recomputed and 
some computations are independent of spin). However, for simplicity we do 
not incorporate this optimisation in our current implementation. 

^http : / /www . fftw.orq/ 

^ ^http : / / www . jasonmcewen .org/ 

" It IS well know that Gauss-Legendre quadrature may be used to construct 
an exact sampling theorem on the sphere. Gauss-Legendre quadrature with P 
points is exact only for a polynomial integrand of order less than or equal 
to 2P — 1. Since neither the associated Legendre functions nor the Wigner 
d-functions are polynomials in cos 8, it does not follow immediately that 
Gauss-Legendre quadrature results in an exact harmonic transform for scalar 
and spin signals on the sphere band-limited at L; nevertheless, this is indeed 
the case. We prove the result for spin s signals on the sphere, thus the scalar 
case will follow simply by setting s = 0. For L samples in 0, we must 
simply prove that the integrand d^^{d) -aGm{9) is polynomial in cosS of 
maximum degree less than or equal 2L — 1. From inspection of |3|, we may 
write the Wigner d-functions as a polynomial of degree £ — s, multiplied by 
(1 - cose)(''-'")/2(l-|- cos 6»)(»+'")/2. Similarly, we may write ^sGmid) 
as a polynomial of degree L — I — s, multiplied by the same factor The 
integrand is therefore polynomial with overall degree L — l + £, which reaches 
a maximum of 2L — 2. Gauss-Legendre quadrature with L samples in 6 may 
thus be used to compute exact spherical harmonic transforms of scalar and 
spin functions band-limited at L. 




Fig. 2. Number of samples A'^ required to represent exactly a signal on 
the sphere of band-limit L for the following sampling theorems: Gauss- 
Legendre sampling theorem (blue/dot-dashed line); DriscoU & Healy sampling 
theorem (green/dashed line); and the sampling theorem developed in this 
ailicle (red/solid line). The inset shows very low band-limits, where the 
difference between Gauss-Legendre sampling and our sampling can have a 
large impact. 



represent a band-limited function, numerical precomputation 
(or lack thereof), the recursions used to compute Wigner 
(i-functions, and numerical accuracy and computation time. 
Finally, we evaluate our sampling theorem in the context of 
potential applications. 

A. Sampling 

The number of samples required to represent a band-limited 
function on the sphere exactly is the fundamental property of 
any sampling theorem, with fewer samples desired. Both the 
Gauss-Legendre and our sampling theorems require in general 
Nqi^ ~ ^MW ~ samples, while the Driscoll & Healy 
sampling theorem requires A^dh ^ 4L^ samples. We therefore 
provide a reduction in number of samples by a factor of two 
compared to the canonical equiangular sampling theorem on 
the sphere. Furthermore, we require A^gl — Nmw = 3(L — 1) 
fewer samples that the Gauss-Legendre sampling theorem, 
which for small band-limits can be significant (as discussed 



in Sec. V-Ei. The optimal number of samples attainable by a 
sampling theorem on the sphere is given by the degrees 
of freedom in harmonic space. No sampling theorem on the 
sphere reaches this bound in general. However, both the Gauss- 
Legendre and our sampling theorem reach the bound in 
the limiting case of small L (we reach the bound for the cases 
L £ {1, 2}, the Gauss-Legendre sampling theorem reaches the 
bound for L = 1 only, while the Driscoll & Healy sampling 
theorem never reaches the bound). In Fig. [2] we plot the 
number of samples against band-limit for various sampling 
theorems. We also plot the position of samples for each of 
these sampling theorems in Fig. [3] 

B. Precomputation 

It is possible to reduce the computational burden of our 
fast algorithms by precomputing the Wigner d-functions for 
argument tt/2 and for all required harmonic indices. Such 
a precomputation would require 0{L^) storage, similar to 
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(a) View of North pole 




(b) View of South pole 



Fig. 3. Sampling schemes for the exact representation of a signal band- 
limited at L = 12. Sample positions are shown for the following sampling 
theorems: Gauss-Legendre sampling theorem (blue dots); DriscoU & Healy 
sampling theorem (green dots); and the sampling theorem developed in 
this article (red dots). Notice that the Driscoll & Healy sampling theorem 
requires approximately twice as many samples on the sphere as the alternative 
samplings. 

the storage requirements of the precomputation for Driscoll 
& Healy based algorithms [26], [27], [29]. For these latter 
approaches precomputation is essential to recover the fast 
algorithms with complexity below 0{L^). However, for high 
band-limits the storage requirements become impractical at 
present (recall that the precomputation at L = 4096 is 
expected to require approximately 77GB of storage). The 
Wigner d-functions may be evaluated accurately and rapidly 
using recursion formulae; therefore to avoid storage problems 
we do not perform any precomputation and instead compute 
Wigner d- functions on-the-fly using the method of Risbo [40]. 

C. Computing Wigner functions 

Since our algorithms require Wigner d-functions evaluated 
on the entire {m,m') plane for argument 7r/2 only, they are 
flexible with regard to the choice of recursion used to compute 
the Wigner plane for a given £. For example, for the on-the-fly 
computation of Wigner d-functions we may use the recursion 
of Risbo [40] (which requires the entire plane) or Trapani & 
Navaza [41] (which is restricted to argument 7r/2), which are 
both of complexity 0{L'^), without altering the overall 0{L^) 
complexity of our algorithms. However, this is not the case for 



alternative algorithms. 

To compute spherical harmonic transforms for the Gauss- 
Legendre and Driscoll & Healy sampling theorems using 
the 0{L^) algorithms described previously, it is necessary 
to compute Wigner d-function values for a single row of 
the Wigner plane only, but for all i and all values of 9. 
For the overall algorithms to remain 0{L'^), the on-the-fly 
computation of the row of the Wigner plane must be performed 
in 0{L) computations. This precludes the use of the recursions 
devised by Risbo [40] and by Trapani & Navaza [41], for 
example, both of which would result in an overall algorithm 
with complexity 0{L'^). Instead, alternative recursions must 
be used, such as the three-term recursion in £ that goes 
pointwise through the Wigner plane (see e.g. (4.5) of [28]; this 
recursion is used in our implementation of the Gauss-Legendre 
sampling theorem). The inflexibility of these algorithms with 
regard to the choice of recursion used to compute Wigner d- 
functions becomes important when we study the stability of 
these recursions in the following section. Of course, this issue 
may be resolved by precomputing Wigner d-functions using 
any recursion, but as we have seen this becomes problematic 
at large band-limits. 

D. Numerical accuracy and computation time 

We evaluate the numerical accuracy and computation time 
of our algorithms that implement our new sampling theo- 
rem, comparing them to our optimised implementation of the 
Gauss-Legendre sampling theorem and to the semi-naive algo- 
rithm [27] in SpharmonicKi tp] implementing the Driscoll 
& Healy sampling theorem. For all cases we do not perform 
any precomputation since this is infeasible for high band-limits 
(recall that the semi-naive algorithm is the fastest algorithm 
implementing the Driscoll & Healy sampling theorem that 
does not require precomputation). In order to assess numerical 
accuracy and computation time we perform the following 
numerical experiment. We generate band-limited test signals 
on the sphere defined by uniformly random spherical harmonic 
coefficients with real and imaginary parts distributed in the 
interval [—1,1]. An inverse transform is performed to synthe- 
sise the test signal on the sphere from its spherical harmonic 
coefficients, followed by a forward transform to recompute 
harmonic coefficients. Numerical accuracy is measured by 
the maximum absolute error between the original spherical 
harmonic coefficients sfem ™d the recomputed values s/^m, 
i.e. e = maxf.m js/fm — s/?m|- Computation time is measured 
by the round-trip computation time taken to perform the 
inverse and forward transform. All numerical experiments are 
performed on a 2.5GHz Intel Pentium dual core processor with 
4GB of RAM and are averaged over five random test signals. 

The maximum absolute error is plotted against band-limit 
in Fig. |4] for different sampling theorems. High numerical 
accuracy is achieved for all sampling theorems at moderate 
band-limits, with errors on the order of the machine precision 
and increasing approximately linearly with band-limit. For 
the Gauss-Legendre sampling theorem we use the gauleg 
routine of Numerical Recipes [43] to compute Gauss-Legendre 

' ^http : //www ■ cs . dartmouth . edu/ -geelong/sphere/ 
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node positions and weights. This method is based on an initial 
approximation for each node position, followed by an iterative 
refinement based on Newton's method, which is likely to 
explain the slightly inferior error performance of the corre- 
sponding algorithms. Furthermore, the algorithms implement- 
ing both the Gauss-Legendre and Driscoll & Healy sampling 
theorems suffer from their lack of flexibility regarding the 
recursion used to compute Wigner d-functions (or associated 
Legendre functions for the spin ,s = case), which necessitates 
the use of the less accurate pointwise three-term recursion 
in I, rather than more accurate alternatives. Consequently, 
the numerical accuracy of our algorithms is superior to the 
optimised implementations of the alternative sampling theo- 
rems. More importantly, however, both the Gauss-Legendre 
sampling theorem and the semi-naive implementation of the 
Driscoll & Healy sampling theorem go unstable between L = 
1024 and L = 2048, due to the instability of the pointwise 
three-term Wigner recursion. As discussed in Sec. |V-C[ these 
algorithms require a recursion with complexity 0{L), in order 
to remain 0{L^) for on-the-fly computation. To resolve the 
instability of these algorithms it would be necessary to use the 
alternative recursion of Risbo [40], making them 0{L^), or to 
perform a precomputation of Wigner d-functions. Neither of 
these solutions are feasible for high band-limits, hence these 
algorithms are restricted to moderate band-limits (unless an 
alternative pointwise recursion can be found that is stable 
to high band-limits). Due to the flexibility of our sampling 
theorem with regard to Wigner recursion, we are able to use 
Risbo's recursion [40], which is stable to at least L = 4096, 
without altering the complexity of our algorithms. Finally, 
note that for the implementations that support spin transforms 
(the Gauss-Legendre and our sampling theorems), the error 
is identical (to statistical noise) for transforms of real and 
complex signals of different spin. 

The computation time for complex signals is plotted against 
band-limit in Fig. |5] for different sampling theorems. For all 
sampling theorems, computation time evolves as 0{L^) as 
predicted. The semi-naive algorithm is slightly faster than our 
algorithm, which is in turn slightly faster than our optimised 
implementation of the Gauss-Legendre sampling theorem. 
Although we plot performance results for our algorithms 
using the recursion of Risbo [40], we also implemented the 
recursion of Trapani & Navaza [41], which we found to be 
approximately 20% faster than Risbo's approach but which 
goes unstable between L — 2048 and L = 4096. Nevertheless, 
the Trapani & Navaza [41] recursion can be used at band- 
limits at or below L = 2048 to provide a considerable speed 
enhancement. In this case, at band-limit L — 1024 we are 
approximately 25% slower than the semi-naive algorithm, 
but twice as fast as the Gauss-Legendre algorithm. However, 
the semi-naive algorithm applies for scalar functions only 
(and requires approximately twice as many samples on the 
sphere), while the alternative sampling theorems also apply 
directly for spin functions on the sphere. Finally, note that for 
the implementations that support spin transforms (the Gauss- 
Legendre and our sampling theorems), computation time is 
identical (to statistical noise) for transforms of signals of 
different spin, as predicted. 




32 64 128 256 512 1024 2048 4096 

L 

Fig. 4. Numerical accuracy of the algorithms implementing the fol- 
lowing sampling theorems: our optimised implementation of the Gauss- 
Legendre sampling theorem (blue/dot-dashed line); the semi-naive algorithm 
in SpharmonicKit implementing the Driscoll & Healy sampling theorem 
(green/dashed line); and our algorithms implementing the sampling theorem 
developed in this article (red/solid line). 0{L) scaling is shown by the 
heavy black/solid line. The algorithms implementing the Gauss-Legendre and 
Driscoll & Healy sampling theorems go unstable between L = 1024 and 
L = 2048, due to the enforced use of the pointwise thi'ee-term Wigner 
recursion. For the Gauss-Legendre and our sampling theorems, which both 
support spin transforms, the maximum absolute error e is averaged over 
complex signals of spin s 6 {0, 2, 10} and a real spin s = signal, with 
one standard deviation error bars shown (in most cases differences are very 
small and error bars cannot be seen easily). Note that for these cases the 
maximum absolute en'or is identical (to statistical noise) for transforms of 
real and complex signals of different spin. 



100000 




Fig. 5. Computation time of the algorithms implementing the fol- 
lowing sampling theorems; our optimised implementation of the Gauss- 
Legendre sampling theorem (blue/dot-dashed line); the semi-naive algorithm 
in SpharmonicKit implementing the Diiscoll & Healy sampling theorem 
(green/dashed line); and our algorithms implementing the sampling theorem 
developed in this article (red/solid line). C'(L'*) scaling is shown by the 
heavy black/solid line. The algorithms implementing the Gauss-Legendre and 
Driscoll & Healy sampling theorems go unstable between L = 1024 and 
L = 2048, due to the enforced use of the pointwise three-term Wigner 
recursion. For the Gauss-Legendre and our sampling theorems, which both 
support spin transforms, the computation time t (seconds) is averaged over 
complex signals of spin s S {0, 2, 10}, with one standard deviation error 
bars shown (in most cases differences are very small and error bars cannot 
be seen easily). Note that for these cases the computation time is identical (to 
statistical noise) for transforms of signals of different spin. 
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E. Applications 

We discuss three potential applications of our sampling the- 
orem in the fields of cosmology, neuroscience and compressive 
sampling (CS). For each case we highlight the enhancements 
that our sampling theorem will afford. 

For the analysis of CMB observations, our sampling theo- 
rem and associated algorithms provide the ability to perform 
fast spherical harmonic transforms that are exact at the very 
high resolution of current and forthcoming CMB observations. 
Furthermore, we can compute harmonic transforms of both the 
temperature and polarisation of the CMB for identical cost. 

The number of samples required to represent a band- 
hmited signal is not only of theoretical interest but also has 
important practical application. A sampling theorem requiring 
fewer samples means a band-limited function can be mea- 
sured exactly for lower cost. This is particularly important 
in appUcations where the cost of acquiring a single sample 
is large. In neuroscience, for example, diffusion magnetic 
resonance imaging (MRI) [44] is one such apphcation, where 
cost is measured in terms of acquisition time. Diffusion 
MRI has received considerable attention recently as a non- 
invasive technique to image structural neuronal connectivity 
in the brain. In this setting, most recent acquisition strategies 
consider sampling on multiple spherical shells for each voxel 
of the brain, from which an orientation distribution function 
(ODF) describing the probability density of neuronal fibre 
directions is recovered [45], [46]. The ODFs of each voxel 
are then combined to recover neuronal connectivity. Given 
the millions of voxels generally considered, at present this 
imaging modality remains too time consuming for clinical use. 
Typically, very low band-limits of order L ~ 10 are considered 
for each spherical shell. Consequently, when adopting an 
exact sampling theorem even the small reduction in number 
of samples between our sampling theorem and the Gauss- 
Legendre approach of A^gl — -^mw = 3(i — 1) can have a 
large impact on the total cost of acquisition. For a typical 
example with acquisitions made on three concentric spherical 
shells of increasing radius, it has been shown that band-limits 
of three, five and nine, respectively, are required to limit 
aliasing to acceptable levels [47]. For such an acquisition, 
the total number of samples, and thus total acquisition time, 
would be reduced by a factor of 13% when replacing Gauss- 
Legendre sampling with our sampUng theorem. This type of 
enhancement is of considerable importance in order to make 
diffusion MRI accessible for clinical use. 

The recently developed theory of CS states that it is possible 
to acquire sparse or compressible signals with fewer samples 
than standard sampling theorems would suggest [48], [49]. In 
these settings, the ratio of the number of required measure- 
ments to the dimensionality of the signal scales linearly with 
its sparsity [48]. By reducing the dimensionality of the signal 
in the spatial domain, our samphng theorem will enhance 
the performance of CS reconstruction on the sphere when 
compared to alternative sampling theorems. Furthermore, for 
sparsity priors defined in the spatial domain, such as signals 
sparse in the magnitude of their gradient, sparsity is also 
directly related to the sampling of the signal. For this class of 



signals, we therefore expect to see an additional enhancement 
in CS reconstruction performance when adopting our sampling 
theorem. The use of CS techniques on the sphere is likely to 
have wide-spread application for a wide range of problems, 
including more efficient acquisition, denoising and deconvo- 
lution on the sphere. In particular, all of these problems are 
faced in diffusion MRI and in analysing the CMB, which we 
are studying currently to evaluate in detail the enhancements 
provided by our sampling theorem. 



VI. Conclusions 

We have developed a novel sampling theorem on the sphere, 
with corresponding fast algorithms, by associating the sphere 
with the torus through a periodic extension. To represent 
a band-limited signal on the sphere exactly our sampling 
theorem requires less than half the number of samples required 
by other equiangular sampling theorems on the sphere and 
an asymptotically identical, but smaller, number of samples 
than the Gauss-Legendre sampling theorem on the sphere. 
The complexity of our algorithms to compute both forward 
and inverse transforms is 0{L^), with identical scaling to a 
standard separation of variables. However, the continual use 
of FFTs reduces the constant prefactor associated with the 
asymptotic scaling considerably, resulting in algorithms that 
may be used to compute harmonic transforms rapidly. Nu- 
merical experiments have shown our algorithms to be approx- 
imately twice as fast as optimised algorithms implementing 
the Gauss-Legendre sampling theorem but approximately 25% 
slower than the semi-naive algorithm, the most universally 
applicable algorithm implementing the equiangular DriscoU & 
Healy sampling theorem. However, the semi-naive algorithm 
applies for scalar functions only, while our sampling theorem 
also applies directly for spin functions on the sphere (whereas 
the computation time for a spin s transform using the semi- 
naive algorithm would scale by 1 + s). Numerical experiments 
have also shown our algorithms to be numerically stable to 
band-limits of L = 4096. Conversely, the algorithms that 
implement the Gauss-Legendre and Driscoll & Healy sampling 
theorems on the sphere are restricted in their use of Wigner 
recursions and, due to the enforced use of the pointwise 
three-term Wigner recursion, go unstable between band-limits 
L = 1024 and L = 2048. 

Our novel sampling theorem and fast algorithms will be of 
practical benefit both at very high and low band-Umits. For the 
analysis of CMB data at very high band-limits (L ~ 4096), our 
sampling theorem yields exact spherical harmonic transforms 
with algorithms that are stable and very accurate. For the 
reconstruction of diffusion MRI images at low band-limits 
(L ^ 10), the reduction in the number of samples required by 
our sampling theorem to represent a band-limited signal may 
be exploited to reduce the cost of acquisition significantly. For 
CS applications in both of these fields and beyond, the reduc- 
tion in dimensionality and sparsity afforded by our sampling 
theorem will enhance the performance of CS reconstruction 
on the sphere. 
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